Федеральная служба государственной статистики ежемесячно оценивает среднюю номинальную заработную плату в России; Известны значения с января 1993 по январь 2016. Необходимо построить прогноз на следующие три года.
Попробуем поделить на число дней в месяце: Ряд не стал более регулярным, так что вернёмся к исходным данным.
STL-декомпозиция ряда:
Оптимальное преобразование Бокса-Кокса и результат его применения:
В данном случае преобразование имеет смысл использовать, так как оно хорошо стабилизирует дисперсию. Попробуем округлить параметр и взять \(\lambda=0\): Результат практически такой же. Далее будем использовать \(\lambda=0\).
## ETS(A,A,A)
##
## Call:
## ets(y = tSeries, lambda = LambdaOpt)
##
## Box-Cox transformation: lambda= 0
##
## Smoothing parameters:
## alpha = 0.8968
## beta = 0.0399
## gamma = 0.1032
##
## Initial states:
## l = 4.6769
## b = 0.0031
## s=0.1716 0.0076 -0.0044 0.0062 -0.0011 0.008
## 0.0304 -0.0293 -0.034 -0.0173 -0.0661 -0.0716
##
## sigma: 0.0335
##
## AIC AICc BIC
## -289.4724 -287.1094 -227.8641
Настроив выбранную модель на обучающей выборке, посчитаем её качество на тестовой:
## ME RMSE MAE MPE MAPE
## Training set 0.02705718 4.218196 2.568493 0.0348505 2.145357
## Test set -20.76670383 30.891186 23.044210 -9.0295990 9.886351
## MASE ACF1 Theil's U
## Training set 0.1821666 -0.06878022 NA
## Test set 1.6343766 0.89928841 1.12668
Остатки:
Достигаемые уровни значимости критерия Льюнга-Бокса для них:
Остатки коррелированы, по всей видимости, модель недостаточно хороша.
Q-Q plot и гистограмма для остатков:
Распределение имеет длинные хвосты и выброс слева.
| Гипотеза | Критерий | Результат проверки | Достигаемый уровень значимости |
|---|---|---|---|
| Нормальность | Шапиро-Уилка | отвергается | 2.120423110^{-16} |
| Несмещённость | Уилкоксона | не отвергается | 0.7450093 |
| Стационарность | KPSS | не отвергается | 0.1 |
Исходный ряд нестационарен (p<0.01, критерий KPSS); сделаем сезонное дифференцирование: Ряд всё ещё нестационарен (p<0.01, критерий KPSS). Проведём ещё одно дифференцирование:
Для полученного ряда гипотеза стационарности не отвергается (p>0.1)
Посмотрим на ACF и PACF полученного продифференцированного ряда:
На ACF значимы лаги 1 и 12, на PACF — 1, 12 и 24. Будем искать модель, оптимальную по AICc, в окрестности ARIMA(1,1,1)(2,1,1)\(_{12}\):
| Модель | AICc |
|---|---|
| ARIMA(1,1,1)(2,1,1)\(_{12}\) | -1044.6245289 |
| ARIMA(0,1,1)(2,1,1)\(_{12}\) | -1040.8500956 |
| ARIMA(1,1,0)(2,1,1)\(_{12}\) | -1042.6126473 |
| ARIMA(1,1,1)(2,1,0)\(_{12}\) | -1027.8636046 |
| ARIMA(1,1,1)(1,1,1)\(_{12}\) | -1046.7032117 |
| ARIMA(1,1,2)(2,1,1)\(_{12}\) | -1042.5910157 |
| ARIMA(1,1,1)(2,1,2)\(_{12}\) | -1043.2935219 |
| ARIMA(1,1,1)(3,1,1)\(_{12}\) | -1045.2735494 |
| ARIMA(0,1,0)(2,1,1)\(_{12}\) | -1033.1606073 |
| ARIMA(1,1,1)(1,1,0)\(_{12}\) | -1015.0624859 |
| ARIMA(1,1,0)(1,1,1)\(_{12}\) | -1044.6899974 |
| ARIMA(0,1,1)(1,1,1)\(_{12}\) | -1042.9255562 |
| ARIMA(1,1,1)(0,1,1)\(_{12}\) | -1048.3684574 |
| ARIMA(2,1,1)(0,1,1)\(_{12}\) | -1043.7285109 |
| ARIMA(1,1,2)(0,1,1)\(_{12}\) | -1042.381999 |
| ARIMA(1,1,1)(0,1,2)\(_{12}\) | -1046.6988068 |
Наилучшая по AICc модель — ARIMA(1,1,1)(0,1,1)\(_{12}\). Посмотрим на её остатки: Видно, что в начале ряда остатки не определены, что логично, поскольку модель сезонная. Отрежем начало ряда остатков и проанализируем их:
Достигаемые уровни значимости критерия Льюнга-Бокса для остатков:
Q-Q plot и гистограмма:
У распределения остатков достаточно тяжёлые хвосты; скорее всего, гипотеза нормальности будет отклонена.
| Гипотеза | Критерий | Результат проверки | Достигаемый уровень значимости |
|---|---|---|---|
| Нормальность | Шапиро-Уилка | отвергается | 2.605680310^{-17} |
| Несмещённость | Уилкоксона | не отвергается | 0.3276633 |
| Стационарность | KPSS | не отвергается | 0.1 |
Из-за отсутствия нормальности в итоговом прогнозе будем использовать предсказательные интервалы, полученные с помощью симуляции.
Настроив выбранную модель на обучающей выборке, посчитаем её качество на тестовой:
## ME RMSE MAE MPE MAPE
## Training set 0.07955178 3.716015 2.369193 0.06590879 2.013870
## Test set -19.42347092 29.977272 21.039086 -8.35565908 9.002061
## MASE ACF1 Theil's U
## Training set 0.1680315 -0.1637931 NA
## Test set 1.4921661 0.9119375 1.096755
Применим функцию auto.arima:
## Series: tSeries
## ARIMA(1,1,1)(2,1,2)[12]
## Box Cox transformation: lambda= 0
##
## Coefficients:
## ar1 ma1 sar1 sar2 sma1 sma2
## 0.7348 -0.5706 -0.6775 0.0262 -0.0534 -0.4508
## s.e. 0.1377 0.1658 0.3990 0.1026 0.3952 0.2579
##
## sigma^2 estimated as 0.001057: log likelihood=528.87
## AIC=-1043.73 AICc=-1043.29 BIC=-1018.7
Предлагается модель ARIMA(1,1,1)(2,1,2)\(_{12}\). Её AICc выше, чем у модели, подобранной вручную. Посмотрим на её остатки:
Отрежем первые 13 отсчётов и продолжим анализ:
| Гипотеза | Критерий | Результат проверки | Достигаемый уровень значимости |
|---|---|---|---|
| Нормальность | Шапиро-Уилка | отвергается | 2.837000110^{-17} |
| Несмещённость | Уилкоксона | не отвергается | 0.389983 |
| Стационарность | KPSS | не отвергается | 0.1 |
Настроив выбранную модель на обучающей выборке, посчитаем её качество на тестовой:
## ME RMSE MAE MPE MAPE
## Training set 0.08046572 3.703396 2.358651 0.06791549 2.004611
## Test set -20.46838493 30.868728 21.933648 -8.78402003 9.374485
## MASE ACF1 Theil's U
## Training set 0.1672838 -0.1581542 NA
## Test set 1.5556116 0.9097194 1.128954
Сравним остатки двух версий аримы, одинаково обрезав их начало так, чтобы у обоих методов они были правильно определены:
##
## Diebold-Mariano Test
##
## data: resres.auto
## DM = 0.36582, Forecast horizon = 1, Loss function power = 2,
## p-value = 0.7148
## alternative hypothesis: two.sided
Критерий Диболда-Мариано не обнаруживает значимого различия между качеством прогнозов.
В целом автоматическая модель сложнее, её остатки не лучше, AICc — выше, а ошибка на тесте больше, так что остановимся на модели, подобранной вручную.
Сравним остатки лучших моделей ARIMA и ETS, одинаково обрезав их начало так, чтобы у обоих методов они были правильно определены:
##
## Diebold-Mariano Test
##
## data: resres.ets
## DM = -2.5201, Forecast horizon = 1, Loss function power = 2,
## p-value = 0.01232
## alternative hypothesis: two.sided
##
## Diebold-Mariano Test
##
## data: resres.ets
## DM = -2.5201, Forecast horizon = 1, Loss function power = 2,
## p-value = 0.006161
## alternative hypothesis: less
Согласно критерию Диболда-Мариано, прогнозы метода ETS значимо менее точные, поэтому в качестве финального выберем прогноз найденной вручную аримы.
Поскольку остатки ненормальны, предсказательные интервалы строятся бутстрепом
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Feb 2016 207.4896 202.2358 213.5791 194.05194 221.2860
## Mar 2016 219.3715 210.2357 230.3308 200.43912 240.0972
## Apr 2016 225.1880 212.1612 240.7841 199.93986 252.4175
## May 2016 223.6507 207.2247 243.2913 193.74074 257.2214
## Jun 2016 232.8283 212.8136 257.3899 195.21792 273.3823
## Jul 2016 224.2908 202.1878 251.8536 182.74209 269.4255
## Aug 2016 216.0811 192.7082 246.1042 170.48145 265.2864
## Sep 2016 218.4019 192.8091 252.1885 166.18473 271.5579
## Oct 2016 220.2080 192.1521 258.0965 161.61769 279.3173
## Nov 2016 220.9240 190.7413 261.2000 158.33086 286.4834
## Dec 2016 284.4404 241.3847 339.8941 201.62677 375.5025
## Jan 2017 208.7604 174.9442 251.5439 143.93969 279.2856
## Feb 2017 207.3798 171.9648 253.8228 140.65478 284.0678
## Mar 2017 219.3534 179.4991 271.5803 145.25588 305.4075
## Apr 2017 225.2421 181.4571 283.3184 146.59201 320.2336
## May 2017 223.7567 178.1807 285.9210 141.75326 324.4977
## Jun 2017 232.9778 183.0249 301.2620 145.22601 344.1172
## Jul 2017 224.4622 174.0226 294.1353 138.00068 338.1499
## Aug 2017 216.2652 166.0780 287.5420 130.64608 328.9515
## Sep 2017 218.6019 166.3323 294.5800 129.86302 339.3105
## Oct 2017 220.4198 164.9945 301.0726 128.74046 348.2100
## Nov 2017 221.1438 163.6116 305.1979 126.97687 355.7924
## Dec 2017 284.7302 207.1066 398.3282 162.21318 465.0166
## Jan 2018 208.9767 150.1895 295.4684 117.20492 345.3028
## Feb 2018 207.5973 148.1065 298.2987 114.57097 350.5245
## Mar 2018 219.5854 155.0142 320.9565 118.87609 378.3988
## Apr 2018 225.4818 157.3582 332.5987 118.83430 397.0380
## May 2018 223.9959 155.1545 336.4681 115.49945 403.6425
## Jun 2018 233.2277 159.0596 354.1100 118.20198 427.0991
## Jul 2018 224.7035 151.5078 346.1786 112.09896 419.3710
## Aug 2018 216.4981 144.2518 337.7193 106.89095 410.0909
## Sep 2018 218.8376 145.3020 346.4828 106.51824 420.5966
## Oct 2018 220.6577 144.6078 352.3042 105.47649 433.9250
## Nov 2018 221.3825 144.2202 357.8302 103.73521 449.3872
## Dec 2018 285.0378 183.8538 466.8318 130.25356 590.6829
## Jan 2019 209.2025 133.5263 348.7313 94.31239 440.9304